arXiv: 1503.0381 Ivl [cond-mat.mes-hall] 12 Mar 2015 


Melting temperature of graphene 


J. H. Los, K.V. Zakharchenko, M. 1. Katsnelson, Annalisa Fasolino* 

Radboud University Nijmegen/Institute for Molecules and Materials, 
Heyendaalseweg 135, 6525A J Nijmegen, The Netherlands and 
* a.fasolino@science.ru.nl 
(Dated: March 13, 2015) 

We present an approach to the melting of graphene based on nucleation theory for a first order 
phase transition from the 2D solid to the 3D liquid via an intermediate quasi-2D liquid. The 
applicability of nucleation theory, supported by the results of systematic atomistic Monte Carlo 
simulations, provides an intrinsic definition of the melting temperature of graphene, Tm, and allows 
us to determine it. We find Tm — 4510 K, about 250 K higher than that of graphite using the same 
interatomic interaction model. The found melting temperature is shown to be in good agreement 
with the asymptotic results of melting simulations for finite disks and ribbons of graphene. Our 
results strongly suggest that graphene is the most refractory of all known materials. 


I. INTRODUCTION 

Surprisingly, understanding of melting is still an open 
problem in condensed matter physics. For instance, af¬ 
ter more than hundred years, there is still no reliable 
theoretical justification of the Lindemann criterion [l|. 
Nevertheless, the melting temperature Tm is well-defined 
as the temperature at which the Gibbs free energy curves 
of the liquid and solid phase intersect. Graphene, how¬ 
ever, is a quasi-2D crystalline membrane 0, which melts 
into a 3D liquid phase Q . The different dimensionality of 
these two phases makes it impossible to determine their 
free energy difference with existing methods, so that it is 
a priori not even clear how to define Tm- Here we show 
that a 2D nucleation theory offers a way to overcome this 
difficulty and give a reliable quantitative value of Tm- 

While for 3D systems, there is no consensus about 
how to describe melting and premelting anomalies, for 
a 2D crystal there is a commonly accepted microscopic 
scenario for melting i: the KTHNY theory [HI. In 
the KTHNY theory, melting occurs via unbinding of 
topological defects, like vortices in superconductors flPj . 
This scenario seems not to apply to graphene. Apart 
from the embedding in 3D space, which affects all its 
structural and thermodynamic properties iP , the na¬ 
ture of defects is completely different from those con¬ 
sidered in the KTHNY theory. In graphene, there are 
neither observations nor predictions of stable single pen¬ 
tagons or heptagon defects (disclinations). Dislocations 
(pentagon-heptagon pairs) have been observed mainly 
at g rain boundaries [l|. The 4.6 eV formation energy 
[ij] of a Stone-Wales (SW) defect made of two adjacent 
pentagon-heptagons (57) pairs is much lower than the 
formation energy of two well separated 57 pairs as oc¬ 
curring in small angle grain boundaries [l3l| . In the Ap¬ 
pendix we show how the formation energy of two 57 pairs 
rapidly increases with separation. The relative stability 
of SW defects make them play a crucial role in the pre¬ 
melting behavior of graphene d and prevents them to 
split into two 57 defects and further, and by this to fol¬ 
low the KTHNY scenario. 

In refd we have shown that, within the same model 


of interactions that we use here, spontaneous melting of 
graphene occurs at a temperature that we call of 4900 
K, giving an upper limit for the true melting temperature 
Tm of graphene. Here we show that the melting of bulk 
graphene can be described by nucleation theory, allow¬ 
ing an unambigous definition of the melting temperature, 
Tm- We find Tm — 4510 K, about 250 K higher than that 
of graphite [l| and sofar the highest of all materials. 

The paper is organized as follows. In section H we 
review previous results d on the spontaneous melting of 
bulk graphene. In section HI we introduce an approach 
based on classical and kinetic nucleation theory to the 
melting in 2D. This approach is applied in section IV to 
bulk graphene. The melting of finite disks and ribbons 
is presented in sections V and VI respectively. Summary 
and conclusions are given in section VH. 


II. SPONTANEOUS MELTING 

In a previous work d: have studied the sponta¬ 

neous melting of graphene by means of Monte Carlo 
(MC) simulations based on the reactive bond order po¬ 
tential LCBOPH M- A suitable Lindemann type order 
parameter for graphene was defined as 712 = (l/o^) < 
{vi — ( 1 / 12 ) >) where a = Ijy/TTp is the atomic 

radius, with p the 2D particle density, is the position 
of the i-th atom and the sum over j runs over the 12 
atoms closest to atom i. For graphene, melting starts 
at 712 ~ 0.1, close to the value found for a strictly 2D, 
triangular lattice [l| . 

In Fig. [1^ we show a snapshot from a MC simulation 
in which a bulk graphene system consisting of N =1008 
atoms with periodic boundary conditions (PBCs) is 
heated at a constant rate. Fig. [Hd and c show the po¬ 
tential energy E and 712 as a function of the applied, 
increasing temperature found in 8 independent MC sim¬ 
ulations. Melting is signalled by a jump in E and 712 , 
but with large variations in the observed melting tem¬ 
perature, Tm, at which this jump occurs for different, 
independent simulations. This jump as well as the co¬ 
existence of solid and liquid parts in Fig. [T^ are typical 
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FIG. 1. a) A snapshot of simulated bulk graphene at the on¬ 
set of melting. The system contains 1008 atoms with PBCs. 
Note the coexistence of solid and liquid parts, b) Energy E 
and c) and order parameter 712 as a function of tempera¬ 
ture T during a linear temperature ramp for 8 independent 
simulations. 


features of first order phase transitions. Obviously, the 
lowest = 4900 K found in Ref. H constitutes only 
an upper limit for the true melting temperature Tm of 
graphene. 


III. NUCLEATION THEORY APPROACH TO 
MELTING IN 2D 

Here we propose a quantitative approach to determine 
Tm by considering the melting as a process initiated and 
dominated by the nucleation of liquid nuclei in the solid 
sheet, a typical scenario for a first order phase transition. 
The possible applicability of nucleation theory is further¬ 
more suggested by our observation that the melting pro¬ 
cess occurs in two steps. First the graphene sheet trans¬ 
forms into a sort of quasi-2D liquid phase consisting of 
entangled and interconnected chains remaining roughly 
within the quasi-2D plane of the previously solid sheet. 
Eventually, the chains start to disentangle and extend to 
3D space, with a diverging simulation box in NPT sim¬ 
ulation at zero pressure applied in this work. In Fig. [T] 
this ’’graphene —>■ quasi-2D liquid —>■ 3D liquid” scenario 
is reflected in the energy curve: first, at melting, the en¬ 
ergy steeply increases by about 1.05 eV after which it 
remains roughly constant for a while before it starts to 
increase further. This two steps scenario via an interme¬ 
diate quasi-2D liquid phase would imply that the melt¬ 
ing temperature of graphene is in fact the temperature at 
which the 2D solid and the quasi-2D liquid phases are in 
equilibrium. Looking at the snapshot in Fig. [T^, it can 
also be argued that before complete melting, when liquid 
nuclei start to form in the solid sheet, the liquid phase 
remains connected with the solid sheet at the edges of the 
nucleus, constraining it to stay within quasi-2D configu¬ 


rations. The importance of assuming the existence of a 
quasi-2D liquid phase with its own thermodynamic prop¬ 
erties is that it allows us to formulate a nucleation theory 
for the melting of a 2D solid embedded in 3D space, the 
case of graphene. If the melting can only take place via 
this intermediate, quasi-2D liquid phase, then this inter¬ 
mediate phase is decisive for the location of the melting 
temperature of graphene. 

For the analysis of melting in terms of nucleation 
theory we performed two series of simulations: simula¬ 
tions at constant temperature and simulations applying 
a linear temperature ramp. To assess the accuracy of 
our approach, we will compare classical nucleation the¬ 
ory (CNT) and kinetic nucleation theory (KNT), as de¬ 
scribed hereafter. We use LCBOPII for the inter¬ 
atomic interactions, as in Ref. Q. We point out that the 
melting temperature of graphite according to LCBOPII 
has been accurately determined to be 4250 K (I^ . 

For the application of CNT to melting in 2D, we write 
the work W (r) done to form a quasi-2D liquid nucleus of 
radius r as: 

W{r) = TTr^pAfj, -b 27rr7s; (1) 

where Ap = pi — ps is the chemical potential differ¬ 
ence per particle between the quasi-2D liquid phase and 
the solid phase, and 7 si is the solid-liquid interface free 
energy. Note that Trr^p is the number of liquid par¬ 
ticles inside the nucleus, which should indeed grow as 
due to the 2D geometry of the solid. For T > Tm, 
Afi < 0 and W (r) has a maximum at the critical nu¬ 
cleus size with radius Tc = —^si/ipAp) equal to Wc = 
W(rc) = —TT^h/^pApi). Using the thermodynamic re¬ 
lation Ap = Ah — T As ~ (Ah/Tm) {Tm — T), with Ah 
and As = Ah/Tm the melting enthalpy (latent heat) and 
melting entropy per particle respectively, the nucleation 
probability is derived to be: 

PnZi' = K exp (-/3IFc) = A exp (2) 

where AT is a kinetic prefactor, /3 = 1 /(/cbT) and we 
defined a = TT^h/{pAh). 

If the melting is dominated by nucleation, then the 
average time t* required for melting to occur is given by 
the solution of the equation: 

Nr^uclin = f P^^J{T{t))dt = 1 (3) 

Jo 

where Nnuci{t) is the number of (super)critical nuclei at 
time t. At constant temperature T > Tm, eq. [ 3 ] simpli¬ 
fies to t* = 1/Pnuci{T). Instead, if the system is slowly 
heated at a heating rate rj, such that T = Tm + pt, we 
have to solve eq. [ 3 ] numerically for t*. The temperature 
at which melting is expected, T^, is then readily obtained 
from Tm=Tm + pt*. 

In principle, CNT does not specify neither the prefac¬ 
tor K nor its temperature dependence. Therefore, we 
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take K constant in our simulation analysis based on 
CNT. Nucleation is then determined by three parame¬ 
ters: a, K and T^. Expressions for the temperature de¬ 
pendence of the kinetic prefactor belong to the domain 
of kinetic nucleation theory (KNT), which allows us to 
perform a more accurate analysis. 

The starting point of KNT 0 is quite different from 
that of CNT. KNT assumes a (liquid) cluster size distri¬ 
bution (CSD) governed by a coupled set of master equa¬ 
tions with appropriate boundary conditions. In CNT, 
the CSD is proportional to exp {—pW{r)) so that it has 
a minimum for clusters of the critical nucleus radius Tc, 
where W is maximal, while it increases exponentially be¬ 
yond the critical size. The latter prediction is clearly 
unphysical. In KNT, this unphysical behavior is avoided 
by requiring the CSD to satisfy the boundary condition 
= 0 , where is the number of clusters that in¬ 
corporate all N atoms of the system. Eventually, KNT 
leads to a stationary state nucleation rate given by: 

rtKNT s rr I P OiTm \ 

(-1) 

where fc is the attachment rate for a particle to join a 
nucleus of the critical size, Z is the so-called Zeldovich 
factor and a as defined before. The attachment rate can 
be expressed as: 

T 

fc = —^^njlexp {-l3Ea) (5) 

J^ref 


where Vref — ksTref/h (with h the Planck constant) 
is an attempt frequency at some reference temperature 
Tref, which we choose to be T^e/ = 4500 K, gy/n^ ac¬ 
counts for the number of perimeter sites around a critical 
nucleus with Uc = {pXgf ) and g a geomet¬ 

rical factor, while Ea is an activation energy barrier for 
the attachment process. For a circular nucleus, the num¬ 
ber of perimeter sites is given by 2Trrcdp = {2dja)y/nf 
where d is the width of one atomic layer (or shell), whence 
g = 2d/a. Following Ref. , the Zeldovich factor Z for 
our 2D case is found to be: 


Z = 


13 

47r aAh ) 


( 6 ) 


Comparing the parameters in KNT to those in CNT, and 
noting that K = Vrefg replaces K in CNT, there are two 
additional parameters: A/i and Ea- However, Ah can 
directly be obtained from the simulations as the melting 
energy (at zero pressure): Ah ~ 1.05 eV (see Fig. [T]d). 
The activation energy barrier Ea is typically a few eV, 
implying that the temperature dependence introduced by 
exp {—(dEa) is rather weak within the temperature range 
of interest and has only a very minor effect (as we ver¬ 
ified). We have taken Ea = 2.6 eV fixed, equal to the 
activation energy for diffusion in a mixed sp-sp^ carbon 
liquid, as found with LCBOPII [l^ . 

LCBOPII include short range (covalent), long range 
(van der Waals) and middle range (MR) interactions [TJ . 



FIG. 2. Comparison of the average, observed melting tem¬ 
perature, Tdn, of bulk graphene (1008 atoms with PBCs) as 
a function of the heating rate obtained from MD (blue dia¬ 
monds) and MC (red triangles) simulations after calibration, 
using LCBOPII without the middle range interactions. Er¬ 
ror bars are based on 48 independent simulations. Lines are 
guides to the eye. 


The latter contribution is a correction that improves the 
reactive properties and, by construction, does not affect 
the equilibrium carbon structures at low T. Since many 
simulations were required to obtain sufficient statistics, 
and the use of the full potential including the MR part 
is computationally demanding for Molecular Dynamics 
(MD) simulations, we have performed MC simulations. 
Following Ref. we then assume that time is propor¬ 
tional to the number of MC cycles Imc (one cycle being 
N trial moves) so that t = CtMC where C is a time cal¬ 
ibration factor assumed to be constant within a limited 
temperature range provided that the acceptance rate is 
kept constant. To assess the validity and accuracy of this 
approximation, we also performed Molecular Dynamics 
simulations, but without the MR part of the LCBOPII 
for the sake of feasibility, and compared the results to 
those from MC simulations using the same potential (i.e. 
LCBOPII without the MR part). The MC and MD sim¬ 
ulations were performed applying a linear temperature 
ramp. In the MD simulations the temperature was raised 
by using the Berendson thermostat and scaling the ve¬ 
locities at every time step, while in the MC simulations 
the temperature was raised every 5 MC cycles. We con¬ 
sidered a series of heating rates, and for each of them the 
average was determined from 48 independent simu¬ 
lations. For each simulation, was determined by the 
intersection of the 712 curve with the horizontal line at 
712 = 0.2 (see Fig. [T^). Plotting the average as a 
function of ry, a very good agreement between the MD 
and MC results was obtained by applying a time calibra¬ 
tion factor equal to C = 5.56 10“"^ ps/cycle, as is shown 
in Fig. [51 We note that for the determination of Tm us¬ 
ing nucleation theory, it is sufficient that t = CtMC holds 
and we do not need C. 
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FIG. 3. Symbols: a) and c) average melting temperature 
as a function of the MC heating rate r]MC and b) and d) MC 
’time’ required for melting at isothermal conditions as 
a function of temperature. The error bars are based on at 
least 48 independent simulations, but those for the isother¬ 
mal simulations at 4750, 4800 and 4850 K are based on 96 
simulations. The solid (red) curves represent simultaneous 
best fits of both sets of simulation data based on predictions 
by CNT (a and b) and KNT (c and d). 




IV. MELTING OF BULK GRAPHENE 

We apply the CNT and KNT described in the previ¬ 
ous section to evaluate of a bulk graphene system 
consisting of N =1008 atoms with PBCs as in Fig. [T] In 
Fig. |3] we present the results of two sets of systematic 
MC simulations, including simulations heating the sys¬ 
tem at a rate rjMC and isothermal simulations. In view 
of the large variations in the respective, observed melting 
temperature and time required for melting, at 

least 48 independent simulations were performed in all 
cases, in order to obtain good statistical averages. Best 
fits to all 13 data points were then determined by using 
the described CNT and KNT approaches, with fitting 
parameters K (or K = Vrejg in the case of KNT), a and 
Tm- In both cases good agreement was obtained, with 
estimates for Tm equal to 4450 and 4510 K according to 
CNT and KNT respectively. This relatively small dif¬ 
ference shows that the exponential factor dominates the 
nucleation rate and that the role of the temperature de¬ 
pendence of the prefactor in KNT is limited. Since KNT 
is more accurate, we conclude that Tm — 4510 K. From 
a we find 7 si = 0.134 eV/A . 


V. MELTING OF FINITE GRAPHENE DISKS 

Searching for verfication of the result from KNT, we 
also performed simulations of the melting of a series of 
roughly circular graphene disks of increasing size. Typ¬ 
ically, their melting temperature, Tm,d, increases as a 
function of cluster size, as was already demonstrated for 
small graphene flakes {N < 55) in Ref. d. According to 
the Pavlov model [2^ specified below, Tjn,d deviates from 
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FIG. 4. a) Snapshot at 4300 K, illustrating edge melting of 
a graphene disk containing 1274 atoms heated at a rate of 
2.5 10“"^ K/cycle. b) Snapshot illustrating the recrystalliza¬ 
tion after cooling the system in a) from 4300 to 3800 K. c) 
Melting curves for different graphene disk sizes: 1) 714, 2) 
1274, 3) 2470, 4) 3820 and 5) 8020 atoms. For the latter clus¬ 
ter, melting curves from 4 independent simulations are drawn 
but they are almost overlapping. The crossing dashed lines 
illustrate the method used to define the melting temperature 
Tm,d of a disk, d) Best fit (red line) according to eq. [7]of the 
average simulated disk melting temperatures (symbols). The 
error bars are based on 4 independent simulations. 


Tm by a finite size correction that vanishes with the in¬ 
verse radius of the disk, allowing for the determination of 
Tm by extrapolation. It should be noticed, however, that 
the analysis of the disk simulations is arguable because 
Tm,d is not rigorously defined and the Pavlov model does 
not consider edge melting. 

For a disk, nucleation is not required and the melting 
starts at the edge, as illustrated in Fig. |3^. According 
to CNT there is no barrier for nucleation at a ID edge, 
which is rough at any T. Therefore there is little variation 
in Tm,d from different simulations, as is illustrated by the 
almost overlapping melting curves from 4 independent 
simulations for the cluster with 8020 atoms in Fig|4j:. 

Interestingly, the melting process is reversible to some 
extent, as long as there are still solid parts in the system. 
In Fig. 0 |d we show the recrystallization that has occurred 
after the system in Fig. 0^ was slowly cooled from 4300 
to 3800 K. 

An appropriate approach to analyse our simulation 
data for disks would be a 2D version of Pavlov’s model 
[13 . Writing the free energy of the solid (s) and quasi-2D 
liquid (/) disk as Gd,p = Npgp + 27ri?p7p„ for P = s,l, 
with Np the number of particles in the solid/liquid disk, 
Rp their radii and 7 p„ the respective edge free energies of 
their interfaces with the vacuum, and imposing equality 
of the chemical potentials through dGd,s/dNs\j^^^j^ = 
dGd,i/dNi\j^^^j^ leads to the following correction for the 
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melting temperature of a finite disk: 


T, 


m,d 


= Trr. 



Ibv - qpiiv ~\ 
PsAhRs ) 


( 7 ) 


where qp = {psj accounts for the 2D density differ¬ 
ence between solid and liquid. 

The results of the disk simulations are summarized in 
Fig. 01 A clear trend in Tm,d as a function of l/i?s is 
observed and a best fit based on eq. [7] yields Tm = 4503 
K, a value very close to that derived from KNT, and 
Isv — qp7iv = 0.58 eV/A . With an (average) value of 
'jsv — 1.05 eV/A for a graphene ribbon from Ref. [ 2 l| 
and taking gp ~ 1 we find jiv ~ 0.47 eV/A . These edge 
energies, and also 7 si = 0.134 eV/A from KNT reported 
above, are physically sound, with 7 s„ > jiy > 7 ^;. 


VI. MELTING OF 2D GRAPHENE RIBBONS 


For further verification of our result for Tm, we have 
also performed constant temperature melting simulations 
of a graphene ribbon. We used a ribbon of a size x Ly 
= 51.12 X 103.23 A^ containing 2016 atoms with PBC 
only in the x-direction. For this geometry, as for the 
disks, the melting proceeds from the free edges, as shown 
in Fig. [^, and there is no nucleation barrier. In this 
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FIG. 5. a) Snapshot of a melting ribbon consisting of 2016 
atom at T = 4500 K. b) The simulated melting speed (sym¬ 
bols), exp {f5Ea)dNi/dt, as a function of the temperature T 
for a graphene ribbon. The inset shows the fraction of liq¬ 
uid particles Ni as a function of Imc for various simulation 
temperatures as indicated. For the completely liquid phase 
Ni ~ 0.5 as many particle still have a 712 < 0.2 due to a 
rather isotropic neighbor surrounding. 

case of a straight edge, the melting speed, dNi/dt with 
Ni the number of liquid particles, which we define as the 
number of particles with 712 > 0 . 2 , is expected to behave 
as K{T)Psi{exp{—PAp) — 1), where Psi — 2Lxd/{'Ka?) 
is the number of sites at the solid-liquid interface and 
k{T) = (i/refT/Tref) exp {-f3Ea). Since |/3A/r| << 1, 
the melting speed can be well approximated as: 

dNi VrefAhexp{-/3Ea)2Lydp ^ ^ 

dt keTrefTm N ^ rn) ( ) 


where we have defined Ni = Ni/N as the fraction 
of liquid particles. In Fig. [5 ]d we have plotted 
exp {pEa)dNi/dt as a function of T obtained from simu¬ 
lations (symbols). The intersection of a best linear fit 
according to Eq. [5] for T > 4500 K with the hori¬ 
zontal axis yields Tm = 4505 K, very close to the pre¬ 
viously found values. It has to be noticed, however, 
that the results for T = 4400, 4450 and 4500 do not 
follow the linear law, as the system actually melts for 
these temperatures. Although finite size effects might 
play a role here, maybe the most important reason is 
that chain-like molecular units evaporate from the liq¬ 
uid edges, as illustrated in Fig. [5^. Since equilibrium in 
the present case is characterized by a liquid edge of fixed 
average width, evaporation must lead to a thinner and 
thinner solid core and finally to complete melting. To 
deal with this issue properly would require to consider 
a three phase equilibrium. From the simulation results 
for T > 4500 K, using Eq. [ 8 ] and the previously found 
time calibration factor t = 5.56 10“'* ps/cycle, one can 
directly determine the attempt frequency Vref- We find 
Vref = 2.3 10** Hz, which is quite reasonable when com¬ 
paring it to ksTref/h = 0.94 10*'* Hz. 


VII. SUMMARY AND CONCLUSIONS 

By analyzing systematic simulations we have shown 
that the melting of pristine bulk graphene, a prototype 
2D solid embedded in 3D space, follows a two-stage sce¬ 
nario of which the first step from the 2D solid to a quasi- 
2D liquid phase is well described by nucleation theory for 
first order phase transitions. As a consequence, graphene 
has an unambiguously defined bulk melting temperature 
Tm, found to be 4510 K using the interatomic potential 
LCBOPH. This value is confirmed by simulations for fi¬ 
nite disks and ribbons, which by extrapolation yield Pm’s 
close to that from nucleation theory, and show reversibil¬ 
ity to some extent. 

Our finding that Tm is higher (about 250 degrees) than 
for graphite is not likely to be an artifact of LCBOPH. As 
a qualitative explanation we suggest that the increased 
stability of graphene is due to a significant positive en¬ 
tropy contribution for the solid phase due to the pos¬ 
sibility of relatively large out-of-plane fluctuations (rip¬ 
pling) at low energy costs, as compared to the situation 
in graphite. 

We expect that the melting scenario found for 
graphene applies to any covalent 2D material. This ex¬ 
pectation is supported by recent melting simulations for 
2D M 0 S 2 , showing coexistence of the 2D solid and liquid 
phase during melting [^ . as we observed for graphene. 
While measuring Tm for graphene seems to be a real chal¬ 
lenge, measurement of Tm becomes easier for 2D material 
with a (much) lower Tm than that of graphene. 
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VIII. APPENDIX: FORMATION ENERGY OF 
STONE-WALES AND 57 DEFECTS. 

A 5775 Stone-Wales (SW) defect can be created in a 
perfect graphene layer by rotating a single bond by 90 
degrees, followed by a geometrical relaxation of the sys¬ 
tem (see Fig. |B^). It consists of two adjacent 57 defects. 
By subsequently rotating another appropriate bond by 
90 degrees (and relaxing the system), the 5775 SW de¬ 
fects can be split into two 57 defects separated by one 


pair of hexagons, as shown in Fig. [BJa. This process can 
be repeated to separate the 57 defects more and more, as 
shown in Fig. [6]: and d. 

The total formation energy as a function of the sepa¬ 
ration between the (57) defects according to LCBOPII is 
given in Fig. [7] and shows that it strongly increases with 
distance. Similar behaviour was previously reported in 
Ref. [l^, but in that case the 57 defects were created at 
grain boundaries in a polycrystalline sheet. In the latter 
work it was shown that the formation energies accord¬ 
ing to LCBOPII are in very good agreement with those 
calculated ab initio within DFT. The strongly increasing 
energy implies that the pathway to disorder by the split¬ 
ting of a SW defects into two 57 defects and their subse¬ 
quent further diffusion, a scenario in agreement with the 
KTHNY theory, is very unfavorable. Indeed, we never 
observed separated 57 defects in our simulations. 
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FIG. 6. Step by step separation of the 5775 SW defect (a) into two 57 defects with one (b), two (c) and nine (d) hexagonon pairs 
in between. Panels e) and f) are side views of the configuration in a) and d), showing much larger out-of-plane deformations 
for the case of separated 57 defects, due to strong stress in that case. 



separation between 57 defects (in number hexagon pairs) 


FIG. 7. Total formation energy as a function of the separation 
between 57 defects according to LGBOPII, calculated for the 
same sample as used for Fig|6] 






